1D Solver

This notebook solves the one dimensional Schrodinger equation for a periodic potential in an interval 0<x<L.  It discretizes the problem on a mesh in real space and then uses the eigensolver in Mathematica.  In this problem the length of a unit cell is taken to be 1 (one), so that the length of the system is equal to the number of unit cells.

Packages:

Needs["ErrorBarPlots`"];

Constants and Definitions:

npercell=10;  (* Number of mesh points in a unit cell *)
ncells=20;   (* Number of unit cells in the lattice *)
nmax=ncells*npercell;    (*  Total number of points in mesh *)
dx=1./npercell;     (* Distance between mesh points *)
invdxsq=1./(dx^2);
V0=128.;    (*  Strength of potential  *) 
p=1.0;
V[x_]:=V0* Cos[2π p x];
Plot[V[x],{x,0,ncells},PlotPoints-> nmax]

Calculation of Eigenvalues

Here we set up the matrix. Note that we enforce periodic boundary conditions.

Clear[hmat];
hmat=Table[0,{i,1,nmax},{j,1,nmax}];
Do[hmat[[i,i]]= 2*invdxsq -V[i*dx],{i,1,nmax}];
Do[hmat[[i,i+1]]=-invdxsq ,{i,1,nmax-1}];
Do[hmat[[i+1,i]]=-invdxsq ,{i,1,nmax-1}];
hmat[[1,nmax]]=-invdxsq;
hmat[[nmax,1]]=-invdxsq;

This cell calculates the eigenvalues and eigenvectors and stores them in a handy format.  The eigenvectors are not sorted originally.

Clear[eigenout,eigval,eigveclist,sortedeigs];
eigenout=Eigensystem[hmat];
eigval=eigenout[[1]];
eigveclist=eigenout[[2]];
sortedeigs=Sort[eigval];
eigvec[n_]:=eigveclist[[n]];

ListPlot[eigval]

ListPlot[sortedeigs]

ListPlot[Take[sortedeigs,21]]

The cell below will plot the the eigenvector associated with the npos eigenvalue, that is if npos=1, it plots the groundstate, if npos=2, it plots the first excited state, etc.

npos=19;
eig=sortedeigs[[npos]];
vecpos=Position[eigval,eig][[1,1]];
vecplot=eigveclist[[vecpos]];
ListPlot[vecplot,Joined->True,PlotRange->All]

ß